      SUBROUTINE  RDCOEF(C,M,N,IPO,AINT,JWD,D,A,PI,SIGMA,DSUM,FCTR)     00010017
CCC                                                                     00020000
CCC           EIGENVALUES AND EXPANSION COEFFICIENTS OF PROLATE &       00030000
CCC           OBLATE SPHEROIDAL WAVE FUNCTIONS WITH REAL 'C'.           00040000
CCC           IPO = 1  -  PROLATE SYSTEM ;   IPO = -1  -  OBLATE SYSTEM 00050000
CCC           AINT :  INITIAL VALUE OF THE EIGENVALUE.                  00060000
CCC           A  :  FINAL EIGENVALUE CORRECTED BY BOUWKAMP,S METHOD     00070000
CCC           D  :  EXPANSION COEFFUCIENTS - FLAMMER,S NORMALIZATION    00080000
CCC           DSUM  :  SUMMATION OF THE D-COEFFICIENTS                  00090000
CCC           SIGMA  :  NORMALIZATION FACTOR OF THE ANGLE FUNCTIONS     00100000
CCC           CALCULATION - REAL DOUBLE PRECISION(REAL*8)               00110000
CCC           C & RESULTS - REAL*8                                      00120000
CCC           EXTERNAL FUNCTIONS : FF, FACTMM                           00130000
CCC           ORIGINALLY PROGRAMMED BY S.ASANO IN 1972 
CCC   -------------------------------------- REVISED IN 1977 ---------- 00130010
      IMPLICIT REAL*8(A-H,O-Z)                                          00140000
      PARAMETER  (JMX=150)                                              00150000
      DIMENSION  DN(JMX),UN(JMX),CN(JMX),ALP(JMX),BET(JMX),GAM(JMX),    00160000
     +           D(JMX)                                                 00170000
C                                                                       00180001
      FLOAT(N)=DFLOAT(N)                                                00190000
      BETA(K,I)=FLOAT(I*(I-1))*FLOAT((K+K+I)*(K+K+I-1))/(FLOAT((2*(K+I) 00200000
     +          -3)*(2*(K+I)-1))*FLOAT((2*(K+I)-1)*(2*(K+I)+1)))        00210005
C      
C     CALL ERRSET(207,256,-1,0,0,208)  
C                                                                       00230001
      FCTR= 1.D0                                                        00231000
      IF(C.GE.17.0.AND.M.GE.17)  FCTR= 1.0D+40                          00232001
          IF(IABS(IPO).EQ.1)   THEN                                     00240001
      C2=C*C*IPO                                                        00250001
      DSUM=0.D0                                                         00260000
      H=C                                                               00270000
      A=AINT                                                            00280000
      L=N-M                                                             00290000
      IFIN0=45 + 2.7*H + L/5 - M/2                                      00300000
       IF(H.LE.2.0) FIN0=20 + L/5 - M/2                                 00300003
       IF(IFIN0.GE.JMX-1) IFIN0= JMX-2                                  00310005
      CRIT=1.0E-06                                                      00320000
      IF(H.GE.3.0) CRIT=5.0E-06                                         00340000
      IF(H.GE.15.0) CRIT=1.0E-05                                        00350000
      IF(H.GE.25.0) CRIT=5.0E-05                                        00351021
      CONVG= 1.D-8                                                      00360003
      IF( DABS(A).LE.1.D0) CONVG= 1.D-08                                00370003
CCCC          PARAMETERS - ALPHA,BETA,GAMMA AND RATIO-N                 00380000
      L0=L-2                                                            00390000
        IF(MOD(L,2).EQ.0)  THEN                                         00400001
      INT=2                                                             00410000
      IFIN=2*(IFIN0/2)                                                  00420000
      ALP0=FLOAT(M*(M+1))+0.5*C2*(1.0-FLOAT(4*M*M-1)/FLOAT((2*M-1)*     00430000
     +     (2*M+3)))                                                    00440005
        ELSE                                                            00450001
      INT=1                                                             00460001
      IFIN=2*(IFIN0/2)+1                                                00470000
        END IF                                                          00480001
 2002 DO 50 JR=INT,IFIN,2                                               00490002
      RM= FLOAT(M+JR)                                                   00500000
      DM= FLOAT(M+M+JR)                                                 00510000
      ALP(JR)=RM*(RM+1.D0)+0.5D0*C2*(1.D0- FLOAT(4*M*M-1)/              00520013
     +       ((2.D0*RM-1.D0)*(2.0*RM+3.0)))                             00530013
      BET(JR)=C2*C2*BETA(M,JR)                                          00540000
      GAM(JR)=C2*DM*(DM-1.D0)/((2.D0*RM-1.D0)*(2.D0*RM+1.D0))           00550013
      IF(JR.LE.2) GO TO 50                                              00560000
      UPLMT=BET(JR)/(A-ALP(JR))                                         00570000
       IF(DABS(UPLMT).LE.CRIT) GO TO 1                                  00580001
   50 CONTINUE                                                          00590000
      JR=IFIN+2                                                         00600000
    1 JW=JR-2                                                           00610000
      CRIT=CRIT/10.                                                     00620000
      IF(JW.LT.10) GO TO 2002                                           00630000
      IF(JW.GE.L) GO TO 300                                             00640000
       GO TO 2002                                                       00650000
CCC           CORRECTION FOR EIGENVALUE  **  BOUWKAMP,S METHOD          00660000
  300 ITERAT=0                                                          00670000
      ISWT=0                                                            00680000
      A0=A                                                              00690000
      LC0=L                                                             00700000
      IS0=2                                                             00710000
      CON=0.8                                                           00720013
      IF(IPO.EQ.-1) CON=1.2                                             00730013
C            ITERATION LOOP                                             00740001
    6 ITERAT=ITERAT+1                                                   00750001
        IF(MOD(L,2).EQ.0)  THEN                                         00760001
      UN(2)=A-ALP0                                                      00770000
      LA=2                                                              00780000
        ELSE                                                            00790001
      UN(3)=A-ALP(1)                                                    00800001
      LA=3                                                              00810000
        END IF                                                          00820001
      JWA=JW-2                                                          00830001
      DN(JW)=BET(JW)/(A-ALP(JW))                                        00840000
      DO 60 IR=LA,JWA,2                                                 00850000
      IRA=JWA+LA-IR                                                     00860000
   60 DN(IRA)=BET(IRA)/(A-ALP(IRA)-DN(IRA+2))                           00870000
      LC=L                                                              00880000
       IF(L0.LT.1) LC=L+2                                               00890001
       IF(H.GE.18.0) LC=L+2                                             00900001
       IF(H.GT.35.0) LC=L+4                                             00910001
       IF(ISWT.EQ.IS0) LC=LC0                                           00920001
      DO 70 KR=LA,LC,2                                                  00930000
   70 UN(KR+2)=A-ALP(KR)-BET(KR)/UN(KR)                                 00940000
      U1=UN(LC)                                                         00950000
      U2=-DN(LC)                                                        00960000
      UPRD=(U1/ DABS(U1))*(U2/ DABS(U2))                                00970000
      IF(UPRD.LE.0.0) GO TO 5050                                        00980000
      ISWT=ISWT+1                                                       00990000
      IF(ISWT.EQ.1.AND.(LC-L).LE.2) LC0=LC+2                            01000000
      IF(ISWT.EQ.1.AND.(LC-L).GE.4) LC0=LC-2                            01010000
      IF(ISWT.EQ.IS0.AND.ITERAT.EQ.IS0) GO TO 5060                      01020000
      GO TO 5050                                                        01030000
 5060 A=A0                                                              01040000
       GO TO 6                                                          01050001
 5050 CONTINUE                                                          01060000
      V2=0.D0                                                           01070000
      P2=1.D0                                                           01080000
      DO 90 IQ=LC,JW,2                                                  01090000
      IF( DABS(P2).LT.1.0D-30) GO TO 90                                 01100001
      P2=P2*DN(IQ)*DN(IQ)/BET(IQ)                                       01110000
      V2=V2+P2                                                          01120000
   90 CONTINUE                                                          01130000
      V1=1.D0                                                           01140000
      IF(LC.EQ.LA) GO TO 5                                              01150000
      P1=1.D0                                                           01160000
      LB=LC-2                                                           01170000
      DO 80 KQ=LA,LB,2                                                  01180000
      KK=LA+LB-KQ                                                       01190000
      P1=P1*BET(KK)/UN(KK)/UN(KK)                                       01200000
      IF( DABS(P1/V1).LT.1.D-20) GO TO 5                                01210013
   80 V1=V1+P1                                                          01220000
    5 DELTA=-(U1+U2)/(V1+V2)                                            01230000
      IF( DABS(DELTA).LT.1.D-08) GO TO 4                                01231011
      CRITA= DABS(DELTA/A)                                              01260000
      IF(CRITA.LE.CON) GO TO 502                                        01270000
       GO TO 6                                                          01280001
  502 A=A+DELTA                                                         01290000
         IF(ITERAT.GE.15) GO TO 4                                       01291013
       IF(CRITA.GT.CONVG) GO TO 6                                       01300001
    4 CONTINUE                                                          01310000
C            ITERATION LOOP END                                         01320001
        IF(MOD(L,2).EQ.0)  THEN                                         01330001
      CN(2)=A-ALP0                                                      01340000
      LA=2                                                              01350000
        ELSE                                                            01360001
      CN(3)=A-ALP(1)                                                    01370001
      LA=3                                                              01380000
        END IF                                                          01390001
      JWA=JW-2                                                          01400001
      LB=LC+2                                                           01410000
      CN(JW)=BET(JW)/(A-ALP(JW))                                        01420000
      DO 3020 JA=LB,JWA,2                                               01430000
      JB=JWA+LB-JA                                                      01440000
 3020 CN(JB)=BET(JB)/(A-ALP(JB)-CN(JB+2))                               01450000
       IF(LA.EQ.LC) GO TO 3050                                          01460001
      DO 3030 JC=LA,LC,2                                                01470000
 3030 CN(JC+2)=A-ALP(JC)-BET(JC)/CN(JC)                                 01480000
 3050 CONTINUE                                                          01490000
CCC           EXPANSION COEFFICIENTS - D                                01500000
      M2=M+M                                                            01510000
        IF(MOD(L,2).EQ.0)  THEN                                         01520001
      LD=4                                                              01530000
      D(2)=CN(2)/GAM(2)                                                 01540000
      R=-1.0D0                                                          01550000
       IF(MOD(L/2,2).EQ.0) R=1.0D0                                      01560001
      ENORM=R/FF(M,N)                                                   01570000
      EVEN= FACTMM(M)                                                   01580020
        ELSE                                                            01590001
      LD=5                                                              01600001
      D(3)=CN(3)/GAM(3)                                                 01610000
      LR=(L-1)/2                                                        01620000
      R=-1.D0                                                           01630000
       IF(MOD(LR,2).EQ.0) R=1.D0                                        01640001
      ONORM=R/FF(M,N)                                                   01650000
      ODD=0.5D0*FACTMM(M+1)                                             01660020
        END IF                                                          01680020
      LE=LD-2                                                           01681020
      DO 15 I=LD,JW,2                                                   01690000
      IN=I                                                              01700000
       IF(DABS(D(I-2)).LT.1.D-74) GO TO 13                              01710025
   15 D(I)=CN(I)/GAM(I)*D(I-2)                                          01720000
      JWD=JW                                                            01730000
      GO TO 14                                                          01740000
   13 JWD=IN                                                            01750000
   14 DO 17 I=LE,JWD,2                                                  01760000
      IRM=I+M                                                           01761022
          CALL  SUBFF(M,IRM,FFD,C0FF)                                   01762022
        IF(MOD(L,2).EQ.0)  THEN                                         01770001
      IA=I/2                                                            01780000
      R=1.D0                                                            01790000
       IF(MOD(IA,2).EQ.1) R=-1.D0                                       01800001
      EVEN=EVEN + C0FF*R*D(I)/FFD                                       01820022
        ELSE                                                            01830001
      IB=(I-1)/2                                                        01840001
      R=1.D0                                                            01850000
       IF(MOD(IB,2).EQ.1) R=-1.D0                                       01860001
      ODD=ODD + C0FF*R*D(I)/FFD                                         01880022
        END IF                                                          01890001
   17 CONTINUE                                                          01900000
        IF(MOD(L,2).EQ.0)  THEN                                         01910001
      D0=ENORM/EVEN                                                     01920000
      D(1)=D0                                                           01930000
      DO 25 J=2,JWD,2                                                   01940000
   25 D(J+1)=D0*D(J)                                                    01950000
        ELSE                                                            01960001
      D(1)=ONORM/ODD                                                    01970001
      DO 35 J=3,JWD,2                                                   01980000
   35 D(J)=D(1)*D(J)                                                    01990000
        END IF                                                          02000001
CC        SUM OF D-COEF. , NORMALIZATION CONST. AND FACTOR PI           02010001
      DSUM=0.D0                                                         02020001
      PI=0.D0                                                           02030000
       JA=JWD                                                           02040004
       IF(MOD(L,2).EQ.0) JA=JWD+1                                       02050001
      DO 21 JB=1,JA,2                                                   02060000
      JJ=JA+1-JB                                                        02070000
       J=JJ-1                                                           02080004
       IF(MOD(L,2).EQ.1) J=JJ                                           02090001
      DSUM=DSUM+D(JJ)                                                   02100000
      Q1= 1.D0                                                          02110019
      Q2= D(JJ)                                                         02111013
        IF(M.GE.1)  THEN                                                02120013
      DO 19 LL=1, M                                                     02130013
      Q1=Q1*DFLOAT(J+LL)                                                02140013
   19 Q2=Q2*DFLOAT(J+LL+M)                                              02141013
        END IF                                                          02150007
      PI=PI+Q1/FCTR*Q2                                                  02160019
   21 CONTINUE                                                          02180000
      SIGMA=0.D0                                                        02190000
      DO 31 JJ=1,JWD,2                                                  02200000
       J=JJ-1                                                           02210004
       IF(MOD(L,2).EQ.1) J=JJ                                           02220001
      Q3= D(JJ)                                                         02231018
      Q4= D(JJ)                                                         02232019
       IF(M.GE.1)  THEN                                                 02240014
      DO 38 II=1, M                                                     02250014
      Q3= Q3*DFLOAT(J+II)                                               02260014
   38 Q4= Q4*DFLOAT(J+II+M)                                             02261014
       END IF                                                           02270004
      CONSIG= 2.D0/DFLOAT(M2+J+J+1)*Q4/FCTR*Q3                          02290019
      SIGMA=SIGMA+CONSIG                                                02300000
      IF(J.GT.L.AND. DABS(CONSIG/SIGMA).LT.1.D-25) GO TO 33             02310008
   31 CONTINUE                                                          02320000
   33 CONTINUE                                                          02330000
       RETURN                                                           02340001
C* 30 WRITE(6,1080)                                                     02350014
C1080 FORMAT(1H0,20X,'*  RDCOEF  -  RETURN WITHOUT CALCULATION  *  EIGEN02360014
C    2VALUE IS NOT DETERMINED ' )                                       02370014
C*     RETURN                                                           02380014
           ELSE                                                         02390001
      WRITE(6,1050)                                                     02400001
 1050 FORMAT(1H0,20X,'*  RDCOEF  -  RETURN WITHOUT CALCULATION  *  IPO  02410000
     1SHOULD BE 1 OR -1' )                                              02420000
       RETURN                                                           02430001
           END IF                                                       02440001
      END                                                               02450000
